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Thermodynamics of mono and di-vacancies in barium titanate 
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The thermodynamic and kinetic properties of mono and di-vacancy defects in cubic (para-electric) barium 
titanate are studied by means of density-functional theory calculations. It is determined which vacancy types 
prevail for given thermodynamic boundary conditions. The calculations confirm the established picture that 
vacancies occur in their nominal charge states almost over the entire band gap. For the dominating range of 
the band gap the di-vacancy binding energies are constant and negative. The system, therefore, strives to 
achieve a state in which under metal-rich (oxygen-rich) conditions all metal (oxygen) vacancies are bound in 
di-vacancy clusters. The migration barriers are calculated for mono- vacancies in different charge states. Since 
oxygen vacancies are found to readily migrate at typical growth temperatures, di-vacancies can be formed 
at ease. The key results of the present study with respect to the thermodynamic behavior of mono and 
di-vacancics influence the initial defect distribution in the ferroelectric phases and therefore the conditions 
for aging. 
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I. INTRODUCTION 

Barium titanate is a prototypical ferroelectric mate- 
rial. At ambient pressure and temperatures above 393 K 
it assumes a para-electric, cubic pcrovskite structure. At 
lower temperatures a tetragonal distortion of the unit cell 
is observed which gives rise to ferro-electricity. At even 
lower temperatures further symmetry reductions lead 
to ferroelectric orthorhombic and rhombohedral phases. 
The sequence of phase transitions is the result of sub- 
tle structural and energetic differences Probably the 
most important technological application of barium ti- 
tanate is in thin-film capacitors^ It also serves as an 
end member in several lead-free ferroelectric alloys^ and 
is used in combination with SrTi03 to obtain tunable 
RF devices^ Because of its technological importance and 
its standing as a prototypical ferroelectric, barium ti- 
tanate has been intensively investigated both experimen- 
tally and theoretically. 

Intrinsic and extrinsic point defects are of special im- 
portance in semi-conductors as they control to a great 
extent the electronic properties of these materials. In 
case of barium titanate, this has motivated a consider- 
able number of studies which used conductivity^^— and 
diffusivit y 13 ' 14 measurements to infer information about 
the point defect equilibria in this material. It is further- 
more assumed that oxygen vacancies play a key role in 
ageing and fatigue of ferroelectric materials by impeding 
domain wall motion or by actin g as local disturbances of 
the polarization (see e.g., Refs. Iliirll9l ). For this reason, 
oxygen vacancies and their associates have been inten- 
sively studied in experimentally (see e. g., R efs. 12014221 ) as 
well as theoretically (see e.g., Refs. [23l-l32l). 
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Experimental investigation of point defect properties, 
however, are usually restricted in that they either provide 
bulk information (e.g., conductivity, diffusivity), which 
corresponds to a macroscopic average over the sample 
and is therefore defect unspecific, or very localized in- 
formation (e.g., electron spin resonance, positron anni- 
hilation spectroscopy), which is, however, restricted to 
certain charge states (unpaired spin states) or types of 
defects (open volumes, vacancies)^ 3 - In general, the cor- 
relation between experimental data and specific defects is 
indirect and therefore potentially ambiguous. Quantum 
mechanical calculations on the other hand are capable 
of providing a very detailed picture of individual defects 
and simultaneously sample a whole variety of both differ- 
ent configurations and charge states. They are therefore 
a very powerful tool for verifying model assumptions and 
for providing fundamental insights into the energetics of 
defects as well as their kinetic and electronic properties. 

The purpose of the present work is twofold. First, we 
aim to verify the defect models which have been employed 
to explain the experimental observations (in particular 
conductivity and diffusivity measurements). Second, we 
explore the formation of di-vacancies which - as argued 
above - is an important ingredient for understanding the 
deterioration of ferroelectric switchability. The calcula- 
tions are carried out for the cubic phase which is the 
most stable phase at high temperatures at which most 
ceramics are processed and at which the initial defect 
distributions arc installed. 



II. METHODOLOGY 

A. Computational setup 

Calculations were carried out within density-func- 
tional theory (DFT) using the Vienna ab-initio simula- 
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tion package (VASP)^r— in combination with the pro- 
jector-augmented wave (PAW) method to represent the 
ionic cores and core electrons . 43 ' 44 In order to find the 
most suitable representation for the exchange-correlation 
(XC) functional preliminary calculations were performed 
within the local-density approximation (LDA) 45 ' 46 as 
well as the generalized-gradient approximation (GGA) in 
the Perdew-Burkc-Ernzerhof (PBE) parameterization4£ 
We furthermore considered the effect of treating the Ti-3s 
and Ti-3p electrons as semi-core states. All calculations 
were carried out with 20 /c-points for Brillouin zone sam- 
pling and a plane wave- cutoff energy of 500 eV. For each 
setup the energy-volume curve of cubic barium titanate 
was computed and the equilibrium lattice constant, cohe- 
sive energy, bulk modulus and its pressure derivative were 
subsequently obtained from a fit to the Birch-Murnaghan 
equation of state. Finally, the full band diagrams were 
calculated at the respective zero pressure lattice constant. 

The results of these preparatory calculations are com- 
piled in Table U in comparison with experimental and 
theoretical data from literature. The influence of the Ti- 
3p and Ti-3s states on the properties included in this 
comparison is very small. The most significant difference 
is the increase of the cohesive energies by about 1% upon 
inclusion of the Ti-3s and Ti-3p electrons in the valence. 
Thus, while for some properties the deep Ti-3s and Ti- 
3p can play a crucial role, in the present context their 
effect is expected to be small. For the sake of computa- 
tional efficiency we did, therefore, not include the Ti-3s 
and Ti-3p electrons in the valence. 

The exchange-correlation functional on the other hand 
has a more pronounced impact. Within the GGA the 
bulk modulus is reasonably well reproduced, but the lat- 
tice constant is overestimated while the opposite applies 
for the LDA results. These findings are consistent with 
the results of previous studies ! 34 ' 48 ' 49 In the present study 
we have decided to employ the LDA. 

In order to determine the properties of the phases in 
Table lU we have employed the same computational set- 
tings as for BaTi03 but varied the number of fc-points 
for each material such as to ensure a convergence of the 
total energy better than 1 meV/unit cell. Energy- volume 
curves were calculated for each of these structures allow- 
ing for full internal relaxation. 



B. Defect calculations 
1. Formation energies 

The formation energy AE D of a defect in charge state 
q is given by^2r— 

AE f D = (E D - E H ) + q(E VBM + E F ) - ^ Afafa (1) 

i 

where Ed is the total energy of the defective system and 
Eh is the total energy of the perfect reference cell. The 
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FIG. 1. Stability diagram for cubic barium titanate as deter- 
mined from density-functional theory calculations. The area 
confined between points A, B, C and D is the chemical sta- 
bility range of BaTi03. The line through points C and D 
corresponds to maximally oxygen-rich conditions and an oxy- 
gen chemical potential of Afi® — OeV. Along lines parallel 
to C-D the oxygen chemical potential is constant. The most 
negative value of A/j.0 = AHf/3 — — 6.61 eV is obtained in 
the upper right corner of the diagram. 



second term describes the dependence on the Fermi level, 
Ep, where -Evbm is the position of the valence band max- 
imum. The variation of the formation energy with the 
chemical potentials of the constituents is given by the 
last term. The difference between the number of atoms 
of type i in the reference cell with respect to the defec- 
tive cell is denoted by An*. The chemical potential fa of 
constituent i can be rewritten as fa = fj,^ lk + A fa where 
l^buik (i eno tes the chemical potential of the standard ref- 
erence state and is equivalent to the cohesive energy per 
atom(see Sect. IIIIBj) . 

For the defect calculations we employed supercells with 
40, 60, 90 and 135 atoms equivalent to 2x2x2 to 3x3x3 
unit cells. The formation energies given in Table IPiTl were 
obtained by extrapolating the data calculated for differ- 
ent supercell sizes to infinite dilution (see below). We 
considered both mono- vacancies (Vq, Vbh, Un) as well 
as nearest-neighbor di- vacancies ( Vsa~ Vq, Vt~ V&) tak- 
ing into account a variety of charge states (see Table ITTT]) . 
Due to their size, it is reasonable to assume that host 
metal ion interstitials have very large formation energies. 
This assumption is supported by analytical potential cal- 
culations which have found Frcnkcl defects to have much 
higher energies than Schottky defects^ 4 - In order to con- 
firm these results we calculated the formation energies of 
all interstitials and anti-sites in the neutral charge state. 
Under metal-rich conditions the formation energies for 
interstitials and anti-sites are at least 3.8 eV larger than 
for the lowest neutral vacancy. Under oxygen-rich con- 
ditions the difference is smaller but even for the most 
favorable case (O on Ti anti-site at point D in Fig. [T]) 
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TABLE I. Bulk properties of cubic barium titanate as obtained from experiment and first-principles calculations. US-PP: 
ultrasoft pseudo-potentials; FP-LAPW: full potential-linearized augmented plane waves; TB-LMTO: tight-binding linear muffin- 
tin orbitals; ASA: atomic sphere approximation; LDA: local-density approximation; GGA: generalized-gradient approximation; 
PBE: Perdew-Burke-Ernzerhof parameterization of the GGA; E c : cohesive energy (eV/f.u.); ag: lattice constant (A); Vq: 
equilibrium volume (A 3 /f.u.); B, B'\ bulk modulus (GPa) and its pressure derivative; E%~ T : direct band gap at P-point (eV); 
Eq~ f : indirect band gap measured between points R and P; m*, m* h \ effective electron (hole) mass at the P-point along (100) 
in units of the electron mass. 
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the energy difference is at least 0.9 cV. In the following 
we therefore neglect interstitials and anti-sites and focus 
exclusively on vacancy defects. It should be pointed out 
that both metal and oxygen interstitials can play impor- 
tant roles in several other oxides with more open lattice 
structures (e.g., TiC>2, ZnO, SnC^). In the case of per- 
ovskitc lattices the interstitial sites are, however, much 
smaller and interstitials lead to significant strains which 
cause large formation energies. 

The Brillouin zone integrations were carried out using 
a non-shifted 2x2x2 T-centered mesh which depending 
on the symmetry of the defect configuration is equivalent 
to 4 to 6 fc-points in the irreducible wedge of the Brillouin 
zone. The plane wave cutoff-energy was set to 500 cV. 

DFT calculations of point defect formation energies 
in semi-conductors are subject to certain shortcomings, 
which in order to obtain physically meaningful results 
must be properly taken into account . 55 ! 56 Although the 
band gap is typically underestimated, energy differences 
within the valence band and conduction bands, respec- 
tively, arc usually rather well described. A simple cor- 
rection of the band structure is, therefore, obtained by 
rigidly shifting the valence band (A£Vb) and the con- 
duction band (AEq B ) with respect to each other. The 
correction energy is then given by 

AE bg = n e AE CB + n h AE VB . (2) 

For a given defect and charge state the number of elec- 
trons in the conduction band n e was determined by inte- 



grating the number of occupied conduction band states. 
The number of holes in the valence band rih was similarly 
obtained as the number of empty valence band states. It 
should be noted that this correction only considers the 
effect of the band gap error on the band energy and as- 
sumes rigid levels. It does not take into account elec- 
tronic relaxations which occur if self-interaction effects 
are properly included. 

Since experimentally the band gap of cubic BaTi03 
cannot be determined for T — > K due to the ferroelectric 
phase transition, we estimated its value by extrapolation 
of the data at higher temperatures which yields E^ pt = 
3.4 eV (compare Fig. 10 in Ref. HI). In order to be able 
to correct for the underestimation of the band gap, we 
furthermore assumed the offset of the calculated band 
structure to be restricted to the conduction band, i.e. 
A£Vb = and AE C b = E c * pt - ££ alc . 

Due to the use of supercells elastic and electrostatic 
interactions between the periodic images of the defects 
need to be taken into account. 55 ' 57 Elastic interactions 
scale inversely with the volume, 0(V~ r ), and therefore, 
have been corrected using finite-size scaling. For charged 
defects electrostatic interactions due to the net charge 
moment are present. The corresponding correction term 
can be expressed in the form of a multipole expansion^ 7 - 
The leading term, which describes monopole-monopolc 
interactions, scales with T/ -1 / 3 and can be calculated ex- 
plicitly. The next term is due to monopole-quadrupole 
interactions. It scales with V~ x and can in principle also 
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TABLE II. Bulk properties of Ba, Ti and O and their com- 
pounds in their resp ective ground-states. Experimental data 
from References I64jro6l . E c : cohesive energy (eV/atom); c/a: 
axial ratio; ro: dimer bond length (A); AHf. enthalpy of 
formation (eV/f.u.); other symbols as in Table[I] 
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TABLE III. Formation energies of mono and di-vacancies 
under the chemical conditions indicated in Fig. [T] Note that 
if AfiBn and A/iTi are given Afj,0 is uniquely determined by 
equation @. The charge state, q, of the defect which deter- 
mines the Fermi level dependence of the formation energies 
via equation ([1} is given in the second column. The num- 
ber of electrons occupying conduction band states (n e ) and 
holes occupying valence band states (nh) are relevant for the 
band gap correction via equation © and are given in the 
third column where positive and negative values indicate nn 
and —7ic, respectively. All energies are given in units of eV. 
The finite-size scaling extrapolation error is given in the last 
column. 
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be evaluated explicitly. Since it displays the same scal- 
ing behavior as the elastic interactions, it is, however, 
conveniently corrected using the same finite-size scaling 
procedure ! 55 ! 57 ' 58 In addition, the latter approach has the 
advantage to avoid ambiguities in the calculation of the 
moments of the net charge distribution. It also allows 
to assess the reliability of the correction by means of the 
error of the linear extrapolation of the data. For strongly 
delocalizcd excess charge distributions higher order terms 
might become important which has, however, not been 
observed in the present calculations. 

Since the magnitude of the electrostatic interactions 
in condensed matter is screened, the monopole-monopolc 
correction term requires knowledge of the static dielectric 
constant. For reasons of consistency the latter should be 
taken from first-principles calculations as well. Using a 
similar approach as in the present work Ghosez et al. 
obtained an average value of e = 57 for barium tita- 
nate^^ which was used in the present work. 



2. Transition energies 

If the formation energies of a given defect in charge 
states qi and qi are known, the thermal (equilibrium) 
transition level can be obtained according to 

e= AE f D ( qi ) - AE f D (q 2 ) 

where AE{j(qi) and AEpfa) denote the formation en- 
ergies at the valence band maximum for charge states q\ 
and g 2j respectively. 

3. Migration energies 

The migration barriers for single vacancies {Vq, Vbh, 
Vri) were calculated employing 40-atom supercells cquiv- 
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alent to 2 x 2 x 2 unit cells. Brillouin zone integra- 
tions were carried out using the same fc-point grids as 
for the formation energy calculations and the plane wave 
cutoff-energy was again 500 eV. In order to obtain the 
saddle points, we applied the climbing image nudged elas- 
tic band (CI-NEB) metho d 60 ! 61 which imposes a minimal 
number of constraints on the transition path. Only jumps 
between nearest neighbor sites were considered. 

Unlike formation energies migration energies are ob- 
tained by calculating differences between configurations 
which are structurally and electronically very similar. As 
discussed in detail in Ref. HH the calculation of migration 
barriers is, therefore, much less sensitive to the errors 
described in the foregoing section. For this reason, com- 
parably small supercells should be sufficient and due to 
the similar electronic structure band gap and potential 
corrections need not be applied. 



III. RESULTS 

A. Band structure 

In agreement with full potential-linearized aug- 
mented plane wave^ and pseudopotential plane wave 
calculations,— our results predict the valence band max- 
imum (VBM) to be located at the R-point. On the other 
hand, in previous tight-binding linear muffin-tin orbitals 
calculations^ within the atomic-sphere approximation 
the VBM was located at the r-point. In all of these calcu- 
lations the difference between the highest occupied levels 
at R and T was calculated to be on the order of 0.1 cV or 
less. The effective masses of electrons and holes at the T- 
point, which can be used for instance for a self-consistent 
determination of the Fermi level and the charge carrier 
concentrations,— are given in Table HI 



B. Chemical potentials and stability diagram 

In order to be able to derive defect formation energies 
from total energy calculations, the thermodynamic reser- 
voirs need to be defined. This requires knowledge of the 
cohesive energies of the constituents in their most sta- 
ble conformation. Furthermore, as the stability range of 
barium titanate is restricted by the formation of compet- 
ing phases such as Ti02 or BaO, the formation energies 
of these compounds have to be computed as well (sec 
Sect. Ill B~Tj) . The results of our calculations of the fully 
relaxed structures are compiled in Table [TT] in compari- 
son with experimental data. The overall agreement with 
the reference data is good. In particular, the calculated 
formation energies compare well with experimental data. 

The range within which the chemical potentials of Ba, 
Ti, and O can vary is restricted by the condition 

A MBa + AyUTi + 3A/X0 = AH f [BaTi0 3 ] , (4) 



which determines the outer triangle in the phase diagram 
shown in Fig. [T] Considering the compounds given in 
Table HI1 the chemical potentials are furthermore subject 
to the following constraints 

A/i Ba + Afi < AH f [B&0] 

A^a + 2A^ < AH f [B&0 2 ] 

A/iTi + 2A/i0 < AH f [Ti0 2 }. 

They confine the stability range of BaTi03 to the gray 
shaded area in Fig. [T] Repeating this analysis with the 
experimental data gives a phase diagram in good agree- 
ment with the calculated one. 



C. Defect formation energies 

The results of our calculations for the formation ener- 
gies are summarized in Table IIIII which shows the forma- 
tion energies of mono and di-vacancies for a Fermi level 
at the valence band maximum (Ep = cV in equation 
([1])) and for four representative combinations of chemi- 
cal potentials, which are indicated in Fig. [T] The vari- 
ation with the Fermi level is shown for the same four 
cases in Fig. [2j Due to the large formation enthalpy 
of BaTiOa the formation energies vary strongly between 
the extremal points of the phase diagram. In the metal- 
rich limit (along A-B) oxygen vacancies prevail. They 
have comparably small formation energies and therefore 
should be abundant defects. In both cases the formation 
energies becomes negative for some Fermi level which de- 
termines the so-called pinning energy, £ P j n . This implies 
that under equilibrium conditions the material cannot 
assume a Fermi level which is closer to the valence band 
maximum than e P i n £1 In the oxygen-rich limit cither bar- 
ium (point C) or titanium vacancies (point D) dominate. 
Fermi level pinning now occurs in the vicinity of the con- 
duction band, which implies that the Fermi level cannot 
be pushed arbitrarily close to the conduction band min- 
imum. 

The equilibrium defect transition levels can be deduced 
from the formation energies using equation They are 
presented in an effective band scheme in Fig. El Since va- 
cancies occur in their nominal charge states (V Ti , V Ba , 
Vq) almost over the entire band gap, only the band edges 
are shown. With regard to the effect of the band gap cor- 
rection term given by equation @, it turns out that if no 
band gap corrections are applied the donor transition lev- 
els (q > 0) are near the calculated conduction band edge, 
whereas they are near the experimental conduction band 
edge if the corrections are included. This consistency in- 
dicates that the application of the band gap correction 
described in Ref. |56| is reasonable in the present case. 

The location of the levels near the band edges is 
in agreement with several models which have been de- 
veloped to reproduce the experimental data (see e.g., 
Refs. It is evident from the extrapolation errors 

in Table IIIII that typically DFT calculations cannot pre- 
dict transition levels with an accuracy better than about 
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FIG. 2. (Color online) Variation of defect formation energies with Fermi level for representative thermodynamic conditions 
indicated in Fig. [T] The numbers indicate the charge states. Parallel lines correspond to identical charge states. The solid and 
dashed lines correspond to mono and di-vacancies, respectively. The arrows indicate the position of the Fermi level pinning 
energy under different conditions. 
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FIG. 3. Transition levels for mono and di-vacancies in 
BaTiOa . Only the band edges are shown. The dashed transi- 
tion levels are positioned inside the valence or the conduction 
bands (indicated by the light gray shaded areas) and are only 
included for illustration. The dark grey shaded areas indicate 
the sum of the extrapolation errors for each transition. 
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FIG. 4. (Color online) Binding energies for Vbs,-V@ and 
Vn- Via di-vacancies as a function of Fermi level. The kinks 
correspond to charge transition points of the isolated defects 
(compare Fig. [2] and Fig. [3]). 



0.1 eV. However, even taking into account this limitation 
the present results show clearly the intrinsic defect levels 
to be very close to the band edges (< 0.4 eV). In par- 
ticular, this applies for the titanium vacancy for which, 
based on more simplistic calculations, the defect level has 
been previously calculated to be located — 1.14eV below 
the conduction band maximum^ A finding which is con- 
futed by the present results. 

The binding energy of a di- vacancy is given by the dif- 
ference between its formation energy and the formation 
energies of the isolated mono- vacancies. It is therefore 
independent of the chemical potentials. Since charge 
transitions occur only within about 0.4 eV of the band 
edges (see Fig. the Fermi level position has very lit- 
tle influence (see Fig. 2]) . Over the widest range of the 



band gap the binding energy of the Vri- Vgj-di- vacancy 
is Ef, = — 1.93eV while the value for the VB a -V0-di- 
vacancy is Ef, = — 0.62cV. Only near the band edges the 
attraction are somewhat reduced. Thus, the association 
of metal and oxygen vacancies is energetically strongly 
favored, especially in the case of the Vri- V0-di- vacancy. 

Note that the formation energy of the VB a ~V0-di- 
vacancy is constant along lines which are parallel to A-D 
in the phase diagram (Fig. [1]) corresponding to the re- 
action BaTiOa + Ve a -V0 — > Ti02- Equivalently, along 
lines parallel to B-D the sum of the formation energies 
of oxygen vacancies and Vri- V0-di- vacancies is constant, 
(BaTi0 3 + V T ~V + V ->• BaO). 
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TABLE IV. Calculated migration energies of mono-vacancies 
in units of eV. The temperature ranges above which the de- 
fects become mobile are given in the last column. The nega- 
tive charge states of the titanium vacancy were not considered 
since already the neutral charge state displays a huge barrier 
and, following the trends for the barium and oxygen vacan- 
cies, the addition of electrons can only be expected to further 
increase this value. 

Defect Charge Barrier Onset of mobility 






5.82 


> 2500 K 


-1 


5.96 


> 2600 K 


-2 


6.00 


> 2600 K 





9.84 


> 4300 K 





1.12 


490-590 K 


+1 


0.97 


420-510 K 


+2 


0.89 


390-480 K 



D. Migration energies 

The calculated migration barriers are compiled in Ta- 
ble IIV1 The smallest migration energies are obtained for 
oxygen vacancies and display a weak charge state de- 
pendence. Experimentally the migration barrier for oxy- 
gen vacancies has been determined to be AH m = 0.91 
fRef. |20|). The calculations compare well with this value. 
In particular, the barrier for the doubly positive charge 
state, in which the oxygen vacancy should occur for a 
Fermi level in the middle of the band gap, is in very 
good agreement with this reference value. We point out 
that the migration barriers for the different charge states 
of the oxygen vacancy are very similar to the values for 
cubic lead titanate calculated by Park within DFT.— 

Using the Einstein relation 6Dt = (r 2 ), one can esti- 
mate the temperature above which a defect becomes mo- 
bile by determining the temperature for which 6Dt ex- 
ceeds V ( r2 )min (compare Ref.l62landl68l). The pre- factor 
for the defect diffusivity can be approximated by the 
lowest optical phonon mode at the T-point which gives 
To « 5 THz (Ref. H). If one assumes a typical isochronal 
annealing time of r = 10 min and a mean defect sep- 
aration between y/{r 2 ) min = 100 nm and 1000 nm, one 
arrives at the values which are given in the last column 
of Table IIVI Obviously the only defects, which are fully 
mobile at typical processing temperatures, are oxygen 
vacancies. In contrast, due to their very large activa- 
tion barriers, the migration of metal vacancies is much 
lower even at temperatures close to the melting point 
(~ 1900 K, Ref. HI). 



IV. DISCUSSION 

Experimentally, at low oxygen partial pressures (n- 
type region) and high temperatures (~ 1300 — 1500 K) the 
dependence between the electrical conductivity and the 
oxygen partial pressure is found to be P^J^ which has 



been assigned to doubly charged oxygen vacancies^— 
In accordance the present calculations predict oxygen 
vacancies to be by far the most important defect un- 
der metal-rich conditions (low oxygen partial pressure) 
and to occur in charge state 2+ almost over the entire 
band gap. At somewhat lower temperatures (< 1300 K) 

a transition to a P$J^ dependence is observed. At least 
two different explanations have been discussed in the 
literature^— Either (1) the charge state of the oxygen 
vacancies changes from 2+ to 1+ or (2) accidental accep- 
tor dopants are present in the material. In order for the 
first explanation to be valid, the 2 + /1+ transition level 
of the oxygen vacancy should be located 1.3 eV below the 
conduction band minimum^ However, since the present 
calculations locate this transition just 0.1 eV below the 
conduction band minimum, they provide support for ac- 
cidental acceptor doping as the cause for the change in 
slope. In fact, a more detailed investigation of the re- 
lation between the conductivity and the oxygen partial 
pressure shows that this mechanism can also explain the 
p 1 ^ dependence observed for higher oxygen partial pres- 
sures (p0 2 > 10 _2 atm, p-type region)^ 

The structural and energetic differences between the 
para-electric, cubic phase and the ferroelectric phases 
play a crucial role in determining the magnitude and 
the temperature dependence of ferro-electricity. Since 
BaTiOa ceramics are typically processed at temperatures 
above the cubic-tetragonal phase transition (T c = 393 K) 
and since the onset of mobility ranges given in Table IIVI 
exceed this temperature, the point defect equilibria which 
are established during cooling should correspond to the 
cubic phase. In particular, initially (prior to ageing, 
see Refs. [ID and Ho ) oxygen vacancies should be rather 
randomly distributed over the symmetrically inequiva- 
lent oxygen lattice sites, even if a clear energetic pref- 
erence exists for a particular lattice site (which is for 
instance the case for "c-site" vacancies in tetragonal lead 
titanat o 23 ' 69 ) . In addition, due to the large energy barri- 
ers for metal vacancy migration the distribution of bar- 
ium and titanium, which is established during growth, is 
expected to be largely maintained if the material under- 
goes ferroelectric phase transitions. 

The binding energies for di- vacancies are negative and 
large which implies a strong chemical driving force for 
their formation. Since the migration barriers for metal 
vacancies are large, they are virtually immobile. In con- 
trast, oxygen vacancies are very mobile at temperatures 
> 500 K. Thus, formation of di- vacancies should occur 
readily at typical growth temperatures by metal vacan- 
cies "capturing" diffusing oxygen vacancies. 

Defect complexes such as di-vacancics or impurity- 
vacancy associates carry a dipole moment (see e.g., 
Ref. [30 for a quantitative calculation). On a cubic lattice 
different orientations of these defect dipoles are energet- 
ically degenerate. In the presence of an electric fiel d 15 i 16 
or for non-cubic lattices^ this degeneracy is, however, 
lifted (i.e., the energies for di- vacancy pairs e.g., ori- 
ented along the [001] and [100] axes differs). As argued 
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above, depending on the barriers, the transition temper- 
ature and the cooling rate defect dipoles might not have 
enough time to achieve the orientation with the lowest 
energy. Since metal vacancies are rather immobile, re- 
orientation of these dipoles is much more likely to occur 
by oxygen vacancy jumps. The gradual reorientation of 
defect dipoles is then determined by the barrier for oxy- 
gen vacancy jumps in the first neighborhood of metal 
vacancies. 



V. CONCLUSIONS 

The thermodynamics and kinetics of vacancy defects 
in paraelectric cubic barium titanate have been stud- 
ied by means of density functional theory calculations. 
First, formation, binding and migration energies were 
derived properly accounting for the shortcomings of the 
DFT method. The binding energies of metal-oxygen 
di-vacancies are found to be negative implying that in 
equilibrium under metal-rich (oxygen-rich) conditions all 
metal (oxygen) vacancies are bound in di-vacancy clus- 
ters. While metal vacancies are practically immobile at 
realistic conditions, oxygen vacancies can readily migrate 
at typical growth temperatures. Di-vacancies can, there- 
fore, form if metal vacancies capture one oxygen vacancy. 

It is furthermore confirmed that mono-vacancies oc- 
cur in their nominal (ionic) charge states (V Ti , V Ba , VA') 
over the widest range of the band gap. The dominant 
charge states for the di-vacancies are (Vti - Vja) and 
(VBa-Vg))*. Intrinsic defect levels are confined to a re- 
gion within ~ 0.4 cV of the band edges. This is partially 
at variance with earlier calculations based on a more sim- 
plistic model, which suggested titanium vacancies to act 
as hole traps ^ 

The temperatures at which defects are immobilized are 
higher than or similar to the cubic-tetragonal phase tran- 
sition temperature. The principal findings of the present 
study with respect to the thermodynamic behavior of 
mono and di-vacancies are therefore also relevant for the 
initial defect distribution in the ferroelectric phases. 
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